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NUMERICAL  INTEGRATION  OF  DIFFERENTIAL  EQUATIONS  BY 
THREE  METHODS  --  EQUATIONS  FOR,  DESCRIPTION  OF 
COMPUTER  ROUTINES,  DISCUSSION  OF  RESULTS 


I.  INTRODUCTION 


Our  work  requires  the  numerical  integration  of  second-order  differential 
equations  to  high  accuracy  (one  part  in  10  ,  say)  relatively  close  to  the  epoch, 

and  good  accuracy  (1  part/lO  )  over  a  long  period.  The  methods  used  must 
also  be  efficient,  lest  many  hours  of  computer  time  be  wasted.  To  meet  this 
need  we  have  dealt  with  three  distinct  methods,  which  the  present  memo  will 
describe  in  various  detail  for  our  own  purposes  and/or  external  utilization. 

The  potential  user  should  be  alerted  that  the  methods  to  be  described  have 
been  adapted  to  meet  the  particular  needs  of  our  program,  and  should  not  ex¬ 
pect  our  specific  realizations  to  be  optimum  for  another  kind  of  problem. 

This  element  of  caution  seems  to  be  lacking  in  the  literature,  as  is  clarity. 

In  fact,  I  developed  our  Adams-Moulton  routine  from  no  more  than  a  pair  of 
predictor-corrector  formulas  given  in  an  unpublished  memo  as  quotients  of 
(finite-difference)  operational  polynomials.  When--  after  all  the  difficulties 
had  been  solved  --  I  found  discussions  of  the  method  in  the  literature,  not 
only  did  these  discourses  fail  to  help  in  refining  our  routine,  but  the  difficul¬ 
ties  we  solved  were  not  treated.  Thus  there  is  in  this  area  ample  opportun¬ 
ity  to  apply  one's  own  insight  and  ingenuity. 

II.  METHOD  OF  ADAMS- MOULTON 


This  method  rests  upon  a  particular  pair  of  quadrature  formulas,  and 


n  +  1 


with  the  aid 


(unlike  some  methods;  e.g.  ,  Runge-Kutta)  steps  from  t  to  t 
of  the  solution  at  previous  steps,  say  y^,  y^  j ,  .  . .  y^  p>  i*  e*  >  met^10(^ 

has  "memory".  Also,  the  predicted  value  yn_j_  j  is  ’’corrected"  by  a  formula 
using  yn  y  hence,  this  is  one  of  a  class  of  "predictor-corrector"  methods, 
a  class  which  is  apparently  superior  to  all  other  methods. 


Let  us  make  the  non-trivial  supposition  that  the  solution  and  its  deriva¬ 
tive  are  known  at  a  sufficient  number  of  previous  points,  and  write 
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Here,  the  V's  are  backward  difference,  with  V°  =  1  by  convention.  We  next 
use  this  estimate  of  Yn+  j  to  obtain  j  from  the  differential  equation.  Then 
we  use  a  second  quadrature  formula  involving  y^^  j  to  correct  yn^_  j: 
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If  necessary,  the  second  step  may  be  repeated  --at  the  cost  of  a  second  com¬ 
putation  of  y^  +  j,  and  perhaps  also  the  difference  table. 

The  explicit  form  of  the  quadrature  formulas  have  been  succinctly  de¬ 
rived  by  Glautz  .  We  give  the  formulas  below  and  leave  the  derivation  for 
the  interested  reader. 


y 


n  +  1 


y  +  h 
7n 


V  +  V2  +  V3  +  .  .  . 


y 


i 

n 


y 


n  +  1 


y  +  h 
7n 


y 


n+  1 


One  of  the  tasks  in  coding  this  method  is  that  of  finding  the  coefficients  a^  and 
b.  up  to  high  degree,  a  painful  task  if  the  polynomial  division  be  done  by  hand. 
We  have  derived  recursive  formulae  to  perform  the  division  which  are  well 
suited  to  machine  computation. 
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where  r^  =  l/i.  And 
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It  has  heretofore  been  common  to  compute  and  print  tables  of  such  coefficients 
to  save  other  workers  the  labor.  However,  a  more  current  assumption  is 
that  everyone  has  a  digital  computer,  for  which  the  above  formulas  are  as 
good  as  numbers. 


It  may  be  noticed  that  the  prediction  and  correction  formulae  have  been 
written  in  a  form  which  appears  to  entail,  computationally,  the  formation  of 
difference  tables  of  y1,  followed  by  multiplication  and  summation.  If  a  sec¬ 
ond  correction  be  required,  then  p  more  multiplications  and  additions  would 
be  necessary.  However,  the  formulae  can  easily  be  rewritten  so  as  to  in¬ 
volve  only  the  y1,  not  its  differences. 
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for  prediction. 


for  correction, 
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The  corresponding  computer  calculations  are  more  efficient  than  from  the  dif¬ 
ference  notation.  However,  this  approach  has  a  defect  which  may  or  may  not 
be  serious  for  a  given  problem;  i.e.  ,  the  behavior  of  the  higher  differences 
is  an  excellent  index  for  the  accuracy  and/or  stability  of  the  integration.  It  is 
common  experience  that  increases  in  the  higher  derivatives  of  the  "exact" 
solution,  or  errors  in  computer  programming,  may  cause  an  increase  in  the 
higher  differences  of  y1,  so  that  without  some  remedy  the  process  becomes 
unstable.  In  our  work,  we  accept  the  liability,  and  simply  stop  an  integration 
when  by  some  simple  test  an  instability  has  obviously  set  in. 

When  used  for  integrating  the  second  order  differential  equations  of  mo¬ 
tion  for  a  planet,  the  Adams-Moulton  class  of  methods  appears  to  have  an  in¬ 
herent  difficulty.  Suppose  that  the  orbit  were  an  ellipse;  then  there  would  be 
four  equations  to  be  integrated  in  this  method: 


y'i  = 


y' 

2 

where  (x^,  y^)  is  the  position  and  (x-,,  y the  velocity.  All  such  singly- 
integrating  methods  require  double  application  as  above  to  progress  from  ac¬ 
celeration  to  velocity  to  coordinate.  But  the  peculiar  form  of  the  equations 
of  motion  for  such  methods  suggests  that  our  computer  integrator  could  eas¬ 
ily  become  a  computer  oscillator.  Indeed,  for  a  given  number  of  differences 
in  the  predictor-corrector  formulae,  use  of  a  step  size  greater  than  a  corre¬ 
sponding  threshold  makes  the  system  oscillate.  We  have  found  that  a  second 
correction  alleviates  the  problem,  at  the  cost  of  a  second  evaluation  of  driving 
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function.  However,  since  the  cure  allows  a  step  size  larger  by  at  least  an 
order-of-magnitude,  the  saving  in  computation  time  is  considerable.  It  is 
interesting  that  in  the  double- correction  mode  there  seems  to  be  no  substan¬ 
tial  difference  between  the  first  and  second  corrected  solutions;  this  appar¬ 
ent  paradox  appears  to  illustrate  the  close  relationship  between  accuracy  and 
stability  in  numerical  integration. 

We  have  programmed  the  Adams- Moulton  method  and  applied  it  with  fine 

th 

results  to  the  N-body  integration  routine.  With  11  differences  its  perform¬ 
ance  as  to  accuracy/computation- time  is  far  superior  to  that  of  the  previous 
method  due  to  Nordsieck  (described  below),  with  a  computer  time  saving  of 
about  an  order-of-magnitude. 

Except  for  the  restrictions  of  fixed  step  size  and  equal- interval  tabular 
output  the  method  has  been  programmed  in  a  flexible,  easily-used  form  for 
general  use. 

There  are  a  number  of  techniques  available  for  starting  the  Adams- 
Moulton  method,  which  requires  the  solution  at  a  sufficient  number  of  points 
at  the  epoch  to  compute  the  difference  tables.  We  simply  use  the  Nordsieck 
routine,  on  which  we  can  rely  for  solutions  near  the  epoch  with  accuracy  to 
14  or  15  decimal  places. 

III.  METHOD  CALLED  "ROYAL  ROAD" 

If  formulas  could  be  found  by  which  one  can  doubly  integrate  in  one 

stride  directly  from  acceleration  to  coordinate,  the  just-discussed  problem 

characteristic  of  single  integration  methods  could  be  avoided.  In  fact  such 

formulas  exist,  and  we  follow  Kopal  in  referring  to  them  as  the  Royal  Road 

nd 

method,  a  description  we  prefer  over  JPL's  name  ("2  sum").  Since  we 

have  not  yet  quite  finished  programming  the  method,  we  shall  limit  the  sec- 

2 

tion  primarily  to  a  listing  of  formulas.  We  follow  Hildebrand  in  their  deri¬ 
vation,  but  carry  the  algebra  (for  finding  coefficients)  in  symbolic  rather 
than  numeric  form,  so  as  to  have  the  computer  find  coefficients  for  high 
order  terms.  For  correction, 
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with  the  b.  the  same  as  in  an  earlier  formula  for  the  Adams- Moulton  method* 
1 

To  carry  out  the  squaring  operation  let  u^  =  =  b^.  Then, 
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For  prediction,  we  operate  on  both  sides  with  the  "advancing1'  operator 
E  =  1/(1  -  V): 
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Since  V  -  2z^  j  +  z^  £,  these  formulae  become 
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where  we  have,  for  clarity,  simply  advanced  the  subscript  of  the  correction 
formula  originally  given.  Of  course,  these  formulae  can  be  expressed  in 
terms  of  the  y1^  and  y'j^_|_  j  instead  of  their  differences,  as  before. 

The  use  of  these  royal  road  formulae  in  PEP  raises  a  series  of  questions 
regarding  computation  of  velocities.  For  our  equations  of  motion  the  driving 
function  y"  is  a  function  of  y'  only  in  the  relativity  term,  which  is  always  nu¬ 
merically  much  smaller  than  other  terms.  In  our  equations  for  partial  deri¬ 
vatives  only  the  b /b  (relativity  factor)  contains  "velocity"  terms,  though  of 
course  these  partials  do  not  need  to  be  as  accurate  as  the  coordinates.  The 
velocities  are  also  used  in  two  other  ways,  both  requiring  considerable  accu¬ 
racy:  (1)  in  processing  "probe"  (e.g.  ,  Mariner  V)  data  the  r  must  be  calcu¬ 
lated  to  better  than  1  mm/sec,  for  which  numerical  differentiation  of  r  may 
not  be  adequate;  and  (2)  in  restarting  numerical  integrations,  where  x,  y,  z, 
x,  y,  z  are  all  needed  to  define  the  osculating  orbit. 

If  the  computation  of  velocities  by  integration  of  y"  could  be  avoided,  a 
significant  portion  of  computer  time  needed  for  PLANET  and  N  BODY  integra¬ 
tions  could  be  saved.  Hence,  we  shall  incorporate  the  option  of  calculating 
velocities  by  numerical  differentiation  with  only  prediction  of  y^  +  since 
there  will  not  be  a  second  correction  of  Yn^. 
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To  derive  the  prediction  formula  we  begin,  as  before,  with  the  correction 
formula  . 
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Now  apply  E  to  both  sides. 
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To  express  these  formulae  in  terms  of  tabular  points,  let 
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with  analogous  equations  for  the  predictor. 
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For  restarting  in  the  absence  of  integrated  velocities  it  will  probably  be 
possible  to  use  velocities  computed  from  the  position  and  perhaps  accelera¬ 
tion  in  the  vicinity  of  the  restarting  time.  We  will  give  below  a  suggested 
formula  for  numerical  differentiation  of  position;  more  elegant  and  accurate 
methods  using  the  acceleration  can  certainly  be  derived,  but  may  not  yet  be 
in  the  literature. 

For  accurate  numerical  differentiation  we  prefer  some  kind  of  central 
difference  formula.  It  is  easy  to  compute  the  mean  of  the  previously  given 
formula  and  the  corresponding  formula  using  forward  differences. 
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We  imagine  the  previous  formula  to  contain  superscripts  b  (for  backward  dif¬ 
ferences),  and  write  analogous  formulas  for  forward  differences: 
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Since  A  yn  =  ^  (-1)  (  k  ^  ^n  +  s  k*  ^  ^°^ows  in  a  not- so- simple  way 
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and  hereby  we  can  obtain  y  We  will  discover  by  experiment  whether  it  is 
accurate  enough. 

Since  the  Royal  Road  method  will  be  stable  without  a  second  correction, 
and  by  using  numerically  differentiated  velocities,  we  expect  a  saving  in 
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computer  time  of  better  than  a  factor  of  2,  an  important  improvement  for  our 
projected  very  long  integrations. 

IV.  METHOD  OF  NORDSIECK 

The  third  method  is  due  to  A.  Nordsieck,  and  is  described  in  detail  in 
Ref.  3.  We  have  programmed  and  thoroughly  checked  out  a  double-precision 
’’package”  based  on  the  reference.  In  this  section  we  comment  on  the  method 
and  upon  some  aspects  of  realizing  package. 

The  Nordsieck  method  is  distinguished  by  two  features:  (1)  it  is  self 
starting,  and  (2)  the  step  size  is  automatically  changed  when  needed.  Since 
the  method  is  one  of  the  class  of  predictor-corrector  methods,  which  involve 
the  "memory"  of  the  solution  for  earlier  values  of  the  independent  variable  (t), 
starting  is  naturally  a  problem  because  in  the  beginning  there  is  nothing  to  re¬ 
member!  The  starting  procedure  in  the  Nordsieck  method  is  so  accurate  and 
easy  to  use  that  the  package  is  very  useful  as  a  starting  routine  for  other 
methods,  and  is  so  used.  The  importance  of  a  good  starting  procedure  should 
be  emphasized;  the  "inherited  error"  due  to  inaccuracy  in  computing  the  first 
few  values  of  a  solution  can  easily  dominate  subsequent  errors. 

The  variable- step- size  feature  performs  as  follows:  the  step  size  (h)  is 
such  that  the  equation  is  integrated  over  a  unit  distance  (from  t  =  tj  to  t  = 
t  j  +  1)  with  an  accuracy  not  worse  than  a  specified  value  (the  program  input 
variable  named  EPS).  Changes  in  h  are  accomplished  efficiently  and  stably, 
so  that  even  poorly  behaved  solutions,  such  as  (continuous)  stepwise-linear 
functions,  are  efficiently  found.  In  the  Planetary  Ephemeris  Program  (PEP) 
work  the  method  proved  to  be  especially  useful  for  integrating  the  equations 
of  motion  of  the  asteroid  Icarus,  which,  because  of  its  high  eccentricity  and 
near- collision  with  the  Earth,  required  the  diminution  of  h  by  about  a  factor 
of  10  at  perihelion  and  at  close  approach. 

Analogous  with  the  utility  of  observing  the  high  order  differences  in  other 
methods,  the  behavior  of  h  provides  an  index  to  the  current  nature  of  the  dif¬ 
ferential  equation.  In  particular,  in  a  number  of  instances  the  onset  of  a  rap¬ 
idly  decreasing  h  has  led  to  the  discovery  of  programming  errors  outside  the 
integrating  package. 
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One  fairly  subtle  point  should  be  mentioned  for  the  benefit  of  those  who 
may  wish  to  use  our  package  for  the  Nordsieck  method.  Ordinarily  no  loss  of 
accuracy  is  entailed  when  h  abruptly  diminishes;  indeed,  it  usually  diminishes 
precisely  to  retain  accuracy.  An  exception  occurs  when  h  is  required  tempo¬ 
rarily  to  decrease  by  more  than  a  factor  of  2  in  order  to  compute  a  solution  at 
a  preset  tabular  point.  We  prevented  this  from  occurring  at  any  tabular  point 
after  the  second  (excluding  the  epoch)  by  allowing  h  to  double  only  at  any  other 
tabular  point,  via  logic  involving  L5.  However,  the  user  must  select  the  ini¬ 
tial  interval  (HC)  and/or  the  tabular  interval  (HMX)  so  that  for  his  problem 
this  circumstance  will  not  occur  at  the  first  or  second  tabular  point. 

Because  of  difficult  algebraic  problems  involved  in  deriving  certain  coef¬ 
ficients,  and  because  of  questions  about  the  stability  of  the  variable  h,  the 

th  th. 

method  has  been  programmed  only  for  the  5  and  6  degrees.  While  satis¬ 
factory  for  general  use,  it  turns  out  that  for  our  work  in  planetary  orbits 
other  methods  of  higher  degree  are  much  more  efficient  --  with  the  probable 
exception  of  Icarus  and  other  bodies  with  similar  orbits.  When  the  higher 
degree  differences  of  a  solution  are  large  for  a  given  h,  the  Nordsieck  method 
responds  by  diminishing  h,  whereas  a  method  which  carries  along  these  higher 
differences  can  condone  the  original  h  without  loss  of  accuracy.  It  may  happen 
that  the  increased  computation  time  for  high-degree  methods  (for  a  given  h) 
may  be  more  than  offset  by  the  larger  allowable  value  of  h,  especially  when 
evaluating  the  right-hand  side  of  the  differential  equation(s)  requires  a  large 
share  of  the  gross  computation  time.  Hence,  we  have  found  that  for  the  inte¬ 
gration  of  the  9-body  problem  the  Nordwieck  method  is  more  consumptive  of 
computer  time  by  a  factor  of  (roughly)  10  in  comparison  with  the  Adams- 
Moulton  method,  for  the  given  required  accuracy. 

Some  additional  comment  is  needed  regarding  the  criteria  for  allowing  h 
to  double.  These  criteria  were  carefully  tailored  to  be  efficient  for  our  kind 
of  celestial  mechanics  integrations,  and  may  or  may  not  be  suitable  for  other 
work.  A  note  of  caution  should  be  offered  to  anyone  who  might  wish  to  change 
them.  We  have  found  that  a  single  step  at  a  doubled  interval  when  non-optimum 
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criteria  were  used  could  produce  a  truncation  error,  propagated  through  sub¬ 
sequent  steps,  which  dominated  all  later  errors.  That  is,  with  no  saving  of  time 
for  the  one  step,  accuracy  was  nevertheless  greatly  diminished.  With  proper 
selection  of  criteria  this  does  not  occur. 
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